A few main packages are dedicated to spatial data handling like import/export, display, geoprocessing, spatial stastistics.
rgdal: interface between R and GDAL (Geospatial Data Abstraction Library) and PROJ4 libraries: raster / vector geospatial data formats and coordinate transformation.
sp: classes and methods for spatial data in R.
rgeos: interface between R and GEOS (Geometry Engine - Open Source) library: area, perimeter, distances, dissolve, buffer, overlap, union, contains…
These packages are still widely used.
sf Website: Simple Features for R
First release: October 20, 2016
sp, rgeos and rgdal functionnalities in one package.
Easier data handling, simpler objects.
Tidy data: compatibility with the pipe synthax and tidyverse operators.
Main author and maintainer: Edzer Pebesma (also sp author)
sf objects data structure:
sfReading layer `martinique' from data source `/home/tim/Documents/prz/flo/data/mtq/martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID): 32620
proj4string: +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs
mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add = TRUE, cex = 1.2, col = "red", pch = 20)Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.000 35297.56 3091.501 12131.617 17136.310
[2,] 35297.557 0.00 38332.602 25518.913 18605.249
[3,] 3091.501 38332.60 0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702 0.000 7177.011
[5,] 17136.310 18605.25 20226.198 7177.011 0.000
Simple union:
mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2, border = "red")Aggregation according to a grouping variable:
library(dplyr)
mtq_u2 <- mtq %>% group_by(STATUT) %>% summarize(P13_POP = sum(P13_POP))
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u2), add = T, lwd = 2, border = "red", col = NA)mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2)
plot(st_geometry(mtq_b), add = T, lwd = 2, border = "red")m <- rbind(c(700015, 1624212), c(700015, 1641586), c(719127, 1641586), c(719127,
1624212), c(700015, 1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border = "red", lwd = 2, add = T)mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col = "red", border = "green", add = T)google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)
mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join = st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col = "lightblue")Several solutions are available:
ggplot2 users can have a look to ggplot2 mapping features (geom_sf) that can mix nicely with ggspatial.tmapcartography is based on base graphics and allow most of basic and advanced cartographic representations. Full disclosure: the speaker is the maintainer of cartography.Here we will focus on cartography and do small examples with ggplot2.
cartographylibrary(sf)
# Import geo layer
sm <- st_read(dsn = "../data/rhone/seine_maritime.geojson", stringsAsFactors = F,
quiet = TRUE)
dep <- st_read(dsn = "../data/rhone/dep.geojson", stringsAsFactors = F, quiet = TRUE)
# change projection (lambert93)
sm <- st_transform(sm, 2154)
dep <- st_transform(dep, 2154)
# Import dataset
csp <- read.csv("../data/rhone/data_seine_maritime.csv")
# merge geolayer and dataset
sm <- merge(sm, csp, by = "INSEE_COM", all.x = TRUE)library(cartography)
plot(st_geometry(sm))
propSymbolsLayer(sm, var = "act")
title("Active Population")# Custom map of active population
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
bb <- st_bbox(sm)
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
plot(st_geometry(sm), col = "cornsilk2", border = NA, lwd = 0.5, add = T)
propSymbolsLayer(sm, var = "act", col = "darkblue", inches = 0.5, border = "white",
lwd = 0.7, symbols = "square", legend.style = "e", legend.pos = "topleft",
legend.title.txt = "Labor Force\n(2014)", legend.values.rnd = 0)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce in Seine-Maritime", sources = "Insee, 2018",
author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE,
scale = NULL, north = FALSE)par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
typoLayer(sm, var = "cat", border = "ivory", lwd = 0.5, legend.values.order = c("agr",
"art", "cad", "int", "emp", "ouv"), col = c("#e3b4a2", "#a2d5d6", "#debbd4",
"#b5dab6", "#afc2e3", "#e9e2c1"), add = TRUE, legend.pos = "n")
legendTypo(title.txt = "Dominant socio-professional Category", col = c("#e3b4a2",
"#a2d5d6", "#debbd4", "#b5dab6", "#afc2e3", "#e9e2c1"), categ = c("Agriculteurs",
"Artisans", "Cadres", "Prof. Inter.", "Employés", "Ouvriers"), nodata = F)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce Distribution in Seine-Maritime", sources = "Insee, 2018",
author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE,
scale = NULL, north = FALSE)sm$pcad <- 100 * sm$cad/sm$act
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(sm, var = "pcad", breaks = bks, col = cols, border = "grey80", lwd = 0.4,
legend.pos = "topleft", legend.title.txt = "Share of managers (%)", add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)library(cartogram)
sm_c1 <- cartogram_cont(sm, "act", 3)
sm$pcad <- 100 * sm$cad/sm$act
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(sm_c1, var = "pcad", breaks = bks, col = cols, border = "grey80",
lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Share of managers (%)",
add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)These maps may allow to hide or, at least, diminish the MAUP.
3904.689 m
grid <- getGridLayer(x = sm, cellsize = 3900 * 3900, type = "regular", var = c("cad",
"act"))
grid$pcad <- 100 * grid$cad/grid$act
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(grid, var = "pcad", breaks = bks, col = cols, border = "grey80",
lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Share of managers\n(en %)",
add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)grid$cad100 <- grid$cad * 100
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
smoothLayer(x = grid, var = "cad100", var2 = "act", typefct = "exponential",
span = 4000, beta = 2, breaks = bks, col = cols, legend.pos = "topleft",
mask = st_buffer(sm, 0), legend.title.txt = "Share of managers* (%)", border = "grey90",
lwd = 0.2, add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)
text(x = 518000, y = 6921123, font = 3, cex = 0.8, labels = "*Potential smoothing\n exponential function\n span = 2.5km, beta = 2")ggplot2library(ggplot2)
ggplot() +
geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
geom_sf(data = sm %>% st_centroid(),
aes(size= cad), colour="#E84923CC", show.legend = 'point') +
scale_size(name = "Active population",
breaks = c(100,1000,5000),
range = c(0,20)) +
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(sm)[c(1,3)],
ylim = st_bbox(sm)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(title = "Active population",
caption = "Insee, 2018\nKim & Tim, 2018")ggplot() +
geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
geom_sf(data = sm, aes(fill = pcad), colour = "grey80") +
scale_fill_gradientn(name = "Share of managers (%)",
colours = carto.pal("green.pal", 3,"wine.pal",3),
values=bks/max(bks))+
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(sm)[c(1,3)],
ylim = st_bbox(sm)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(title = "Managers",
caption = "Insee, 2018\nKim & Tim, 2018")CRAN Task View: Analysis of Spatial Data
See file data_prep.R for data extraction.
# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)# Load data
adm <- readRDS("../data/iris_p13.rds")
feat_sir <- readRDS("../data/sir_p13.rds")
feat_osm <- readRDS("../data/osm_p13.rds")
# Define margins
par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
# Plot the maps
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "SIR", sources = paste0(nrow(feat_sir), " restaurants"),
author = "Kim & Tim, 2018", scale = NULL, frame = FALSE, tabtitle = TRUE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, bty = "n",
cex = 0.7, col = "red")
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "OSM", sources = paste0(nrow(feat_osm), " restaurants"),
author = "Kim & Tim, 2018", scale = 0.5, north = T, frame = FALSE, tabtitle = TRUE)Explore datasets interactively
The default is not bad: